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Abstract 

Temporal lobe epilepsy is one of the most common chronic neuro- 
logical disorder characterized by the occurrence of spontaneous re- 
current seizures which can be observed at the level of populations 
through electroencephalogram (EEG) recordings. This paper sum- 
marizes some preliminary works aimed to understand from a theoret- 
ical viewpoint the occurrence of this type of seizures and the origin 
of the oscillatory activity in some classical cortical column models. 
We relate these rhythmic activities to the structure of the set of pe- 
riodic orbits in the models, and therefore to their bifurcations. We 
will be mainly interested Jansen and Rit model, and study the codi- 
mension one, two and a codimension three bifurcations of equilibria 
and cycles of this model. We can therefore understand the effect of 
the different biological parameters of the system of the apparition of 
epileptiform activity and observe the emergence of alpha, delta and 
theta sleep waves in a certain range of parameter. We then present a 
very quick study of Wendling and Chauvel's model which takes into 
account GABAa inhibitory postsynaptic currents. 

1 Introduction 

Epilepsy is a common chronic neurological disorder affecting 1% of the world 
population. It is characterized by the recurrence of seizures that strongly alter 
the patient quality of life. Epileptic seizures are transient manifestations of 
abnormal brain activity characterized by excessive and highly synchronous firing 
in networks of neurons distributing over focal or more extended cerebral regions. 
These networks are often referred to as "epileptogenic networks" (REF). 
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It can manifest as an alteration in mental state, tonic or clonic movements, 
convulsions, and various other psychic symptoms. Different types of seizures are 
distinguished according to whether the source of the seizure within the brain 
is localized (partial or focal onset seizures) or distributed (generalized seizures) 
further divided according to the possible loss of conciousness and to the possible 
effect on the body. 

Indeed, although the exact mechanisms leading to the various forms of 
epilepsy are still largely unknown, it is now commonly accepted that seizures 
result from a change in the excitability of single neurons or populations of neu- 
rons (see [3T]). Indeed, it is known, from epilepsy research, that the nature of 
the interactions between neurons and the properties of neurons themselves are 
altered in epileptogenic networks, although these alterations are not completely 
described and understood. Most of the studies suggest that both hyperex- 
citability and hypersynchronization processes occur in epileptogenic networks. 
However, although these two concepts are now widely accepted, most neurosci- 
entists agree on the fact that they must be better defined and that associated 
mechanisms still have to be quantitatively characterized. Indeed, epileptogenic 
networks are complex dynamical systems in which increased excitability and 
synchronization may be caused by many different nonlinear factors, ranging 
from subcellular to neuronal population level (REF REF). 

In this context, the objective of this paper is show how theoretical stud- 
ies performed in physiologically-plausible computational models of neuronal as- 
semblies ("neural mass models") can allow us to establish some relationships 
between excitability-related parameters in models and some characteristic elec- 
trophysiological patterns typically observed in local field potentials (LFPs) or 
in the EEG recorded under normal or epileptic conditions. 

More particularly, we report results from a mathematical study (stability, 
bifurcation) of two well-established biologically-inspired macroscopic models of 
EEG generation. The first model, refereed to as the Jansen and Rits model 
[291 128] . is a minimal model of a neuronal population comprising two sub- 
populations of cells: pyramidal neurons and interneurons. Despite its relative 
simplicity, this model was shown to simulate signals with realistic temporal dy- 
namics as encountered is real EEG (background activity, alpha activity, sporadic 
or rhythmic epileptic spikes) . The second model, refereed to as the Wendling 
and Chauvels model [49J, [4Sl [51] , is an extended version of the former. It was 
intended to reproduce some architectonic features of the hippocampus, a sub- 
cortical structure very often involved in temporal lobe epilepsy (TLE). This 
model which accounts for three sub-populations of cells (pyramidal neurons and 
two-types of interneurons) was shown to simulate fast oscillations (beta, low 
gamma frequency band) as encountered at the onset of seizures in TLE. This 
fast onset activity was not represented by the Jansen and Rits model. 

Seizures are characterized by single electrical transients called spikes at a 
slow time resolution (hundreds of milliseconds to seconds), macroscopic events 
which have to be distinguished from spikes of single nerve cells, that last only 
1 or 2ms. The EEG represents a set of potential recorded by multiple electrodes 
on the surface of the scalp. The electrical signal recorded is a measure of the ex- 
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tracellular current flow from the summated activity of many neurons, distorted 
by the filtering, attenuated by layers of tissues and bone, representing the ac- 
tivity of neural population of the brain. EEG patterns are characterized by the 
frequency and the amplitude of electrical activity in the range of l?30Hz (some- 
times larger) with amplitudes of 20 ? 100?V, and the frequencies observed are 
divided into four groups: delta (0.5-4 Hz), theta (4-7Hz), alpha (8-13 Hz) and 
beta (13-30 Hz). As neuronal aggregates become synchronized, the amplitude 
of the summated current becomes larger. 



2 Neural mass models 

Jansen's neural mass model was first introduced by Lopes Da Silva and col- 
leagues in 1974, and studied further by Van Rotterdam et al. in 1982 [541 l3"5l . 
These authors developed a biologically inspired mathematical framework to sim- 
ulate spontaneous electrical activities of neuronal assemblies measured for in- 
stance by EEG, with a particular interest for alpha activity. In their model, three 
neuronal populations interact with both excitatory and inhibitory connections. 
Jansen et al. [29] [28] discovered that besides alpha activity, this model was 
also able to simulate evoked potentials, i.e. EEG-like activity observed after a 
sensory stimulation. More recently, Wendling and colleagues used this model 
to synthesize activities very similar to those observed in epileptic patients [?§] , 
and David and Friston studied connectivity between cortical areas with a sim- 
ilar framework [HE]. Nevertheless, one of the main issue of Jansen's model is 
that it is not able to produce all the rhythms in play in epileptic activity or ob- 
served in EEG recordings. This limitation of the model lead Wendling, Chauvcl 
and their colleagues [5T] [3D] to develop an extended Jansen model to better 
reproduce hippocampus activity. This model is based on neurological studies 
[30] revealing that gamma-frequency oscillations were linked with the inhibitory 
intcrneurons in the hippocampal networks, and that two types of GABAa in- 
hibitory postsynaptic currents may play a crucial role in the formation of the 
nested theta/gamma rhythms in the hippocampal pyramidal cells. For these 
reasons, we will also study the 5-populations extended Wendling and Chauvel's 
model. 



2.1 Jansen and Rit's model 



Description of the model The initial Jansen and Rit's model features a 



population of pyramidal neurons (central part of figure 1(a)) that receive exci- 
tatory and inhibitory feedback from local inter-neurons and an excitatory input 
from neighboring cortical units and sub-cortical structures like the thalamus. 
Actually the excitatory feedback must be considered as coming from both local 
pyramidal neurons and genuine excitatory interneurons like spiny stellate cells. 



Figure 1(b) is a block diagram representation of figure 1(a) representing the 
mathematical operations performed inside such a cortical unit. The excitatory 
input is represented by an arbitrary average firing rate p(t) which can be random 



4 




Figure 1: (1(a) (Neural mass model of a cortical unit: a pyramidal cells popula- 
tion interacts with an excitatory and an inhibitory population of interneurons. 
(1(b) I Block representation of a unit, h boxes are synaptic transformations, 
S boxes simulate cell bodies of neurons and is the sigmoidal transform of the 
membrane potential into an output firing rate. The constants Ji account for 
the strength of the synaptic connections between populations. 



(accounting for a non specific background activity) or deterministic, accounting 
for some specific activity in other cortical units. 



The postsynaptic system The functions h e (t) and hi{t) of figure 1(b) are 
the average EPSP and IPSP. They convert the average input firing rate into an 
average excitatory or inhibitory post-synaptic potential. Following the works of 
van Rotterdam et al [13] these transfer function are of type: 



h(t) 



al3te-( 5t t > 
t < 



In other words, if x(t) is the input to the system, its output y(t) is the convo- 
lution product h-kx(t). The parameters a determines the maximal amplitude 
of the post-synaptic potentials and (3 corresponds to the characteristic time of 
integration, which is mainly linked with the kinetics of synaptic transmission 
and with the averaged distributed delays in the dendritic tree. 
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These transfer functions are solutions of the second-order differential equa- 
tion: 

which can be conveniently rewritten as a system of two first-order equations 

z{t) = a(3x(t) - 2az(t) - a 2 y(t) 

These two constants are different for the excitatory and inhibitory popula- 
tions to fit the experimentally recorded EPSP and IPSP functions. The parame- 
ters a and (3 have been adjusted by van Rotterdam |46| to reproduce some basic 
properties of real post-synaptic potentials. Following their works, parameters 
can be set as follows: 



Excitatory population a =: A = 3.25 mV /3=:a = 100s 1 
Inhibitory population a =: B = 22 mV (3 =: b = 50 s~ 1 

Firing rates The activity of the population is considered to be the related 
mean firing rate. This mean firing rate is modelled as a sigmoidal transformation 
of the average membrane potential (see, e.g. [32]). The function Sigm we chose 
in this model approximates the functions proposed by Freeman [15], and has 
the form: 

bigmlv) = — - — (1 + tanh -{v — vn)) = -. r . 

° v i 2 2 1 + e r ( v °~ v > 

where v max is the maximum firing rate of the families of neurons, vq is the value 
of the potential for which a 50% firing rate is achieved and r is the slope of the 
sigmoid at vq. 

The excitability of cortical neurons can vary as a function of the action of 
several substances and vo could potentially take different values. In our model, 
we nevertheless consider a fixed value vq = 6mV as suggested by Jansen on the 
basis of experimental studies due to Freeman [T7] . The works of the latter also 
suggest that v max = 5s _1 and r = 0.56 mV' 1 , the values used by Jansen and 
Rit. 



Interconnections The three neural populations defined interact through ex- 
citatory and inhibitory synapses. The number of synapses established between 



two neuronal populations is denoted by Ji for i = 1 ... 4 as in diagram 1 (b) 



They are considered to be constant and proportional to the average number of 
synapses J between populations. We denote by on the related coefficient (i.e. 
Ji = 014 J). These coefficients can be seen as average probabilities of establishing 
connections between two populations. On the basis of several neuroanatomical 
studies (see for instance [1]) these quantities have been estimated by counting 
synapses, and their numerical values is given in table [TJ 



G 



Parameter 


Interpretation 


Value 


A 


A vpraffp P'veitfltnrv ^vnflrttif (rain 


3 2hinV 


B 


Average inhibitory synaptic gain 


22mV 


1/a 


Time constant of excitatory PSP 


10ms 


1/6 


Time constant of inhibitory PSP 


20ms 


ai, a 2 


Average probability of synaptic contacts in the 


a 4 = 1, ai = 0.8 




feedback excitatory loop 




a 3 , a 4 


Average probability of synaptic contacts in the 


a 3 = a 4 = 0.25 




slow feedback inhibitory loop 




J 


Average number of synapses between populations 


135 


VQ, V m ax, r 


Parameters of the Sigmoid S 


v = 6mV, v max = 5s" 1 






r = OMmV- 1 



Table 1: Numerical values used in Jansen's original model 



Note that we consider here constant synaptic weights. The variability in the 
connectivity weights can be taken into account (see e.g. [11]), and the resulting 
equation is way more complex than the equation we deal with in the present 
article. 



Equations of the model 



the three variables tjq, yi and 1/2 of figure 1(b) 



Following Janscn and Rit's initial work, we consider 
To write the system into a set 



of first-order ordinary differential equation we introduce the derivatives of these 
variables, yo, y\, j/2, noted j/3, j/4 and y$, respectively. We therefore obtain a 
system of 6 first-order differential equations that describes Jansen's neural mass 
model: 



Mt) = Ite(*) lfe(t) = AxSigm[ yi (t) - y 2 (t)} - 2ay 3 (t) - a 2 y (t) 
y\{t) = y 4 {t) y 4 (t) = Aa{p(t) + J 2 Sigm[J iyo {t)}} - 2ay 4 (t) - a 2 yi {t) 
lh(t) = Vs(t) = BbJ 4 Sigm[J 3 y {t)} - 2by 5 (t) - b 2 y 2 (t). 

(1) 

Our study will focus on the variable y = y% — y 2 which models the membrane 
potential of the pyramidal cells since their electrical activity corresponds to the 
EEG signal: pyramidal neurons throw their apical dendrites to the superficial 
layers of the cortex where the post-synaptic potentials are summed, and there- 
fore account for the essential part of the EEG activity [31] ■ Table [T] summarizes 
the different numerical values of the original Jansen and Rit's model. 

Dimensionless reduction of Jansen's model We now change variables in 
order to reduce the number of parameters in the system. First of all, we chose 
the time scale of the excitatory population as our new unit time and define 
t = at. Therefore the new characteristic scale of the inhibitory variable will be 
given by the ratio d = -, G the ratio of postsynaptic amplitudes B/A. The 



7 



dimensionless sigmoidal transform we define is 



where fco = e rv °. We eventually introduce the dimensionless scaled input pa- 
rameter P = (rA) 2- and and the dimensionless scaled connectivity strength 
j = {rA) !W j. The new state variables (Y , X, Y 1 ,Y 2 ,Y 3 ,Y 4 ,Y 5 ) defined by 

Y (t) = Jry Q (T/a) 
Y l {T)=ry l { T /a) i = 1, 2 
X = Y 1 -Y 2 . 

satisfy the equations 



Y 


= Y 3 




X 


= Y 4 -Y 5 




Y 2 


= Y 5 




< . 


= jS(X)-2Y 3 -Y 




% 


= P + a 2 jS(a 1 Y )-2Y 4 - 


(Y 2 + X) 




= d ai GjS(a 3 Y )-2dY 5 - 


~d 2 Y 2 



(2) 



They depend upon the nine dimensionless parameters ctj, i = 1, • • • , 4, fco, d, G, 
j, and P as opposed to the 11 parameters of the original model (see tabled]). 
The numerical values these parameters corresponding to table [T] are given by: 



G 


= | = 6.7692 




<1 


= £ = 0.5 




ai 


= 1 a 2 


= 0.8 


a 3 


= 0.25 a 4 


= 0.25 


log(/c ) 


= rvQ = 3.36 




j 


= (rA) Vm «°> J = 12.285 





(3) 



2.2 Wendling and Chauvel's model of hyppocampus ac- 
tivity 

One of the main drawbacks of Jansen's model is that it is unable to generate 
certain types of cortical activity, for instance seen in epileptic activity. As an 
example, it cannot reproduce the type of fast activity observed at the onset of 
seizures in limbic structures. In order to propose a better cortical mass model, 
Wendling, Chauvel and colleagues [5TJ [SD] revisited the organization of subsets 
of neurons and interneurons, focusing on the hippocampus activity. Based on 
these considerations, they proposed a new neural mass model whose parameters 
were estimated using real EEG signals. Their model is shown as a block diagram 
in figure Fig. [2j The main difference with Jansen's initial model is the addition 
of somatic-projecting inhibitory neurons (GAB Aa, fast receptors). 



8 



External inputs 



Main population 




nterneurons 



(a) Populations and connectivities 



yo 



(t) 



(t) 



yi 



(b) Block diagram 



J\hR t) 



Figure 2: Neuronal population model based on the cellular organization of the 



hippocampus. 2(a) Schematic representation of the model. The pyramidal 
cells population projects to and receives feedback excitatory and inhibitory in- 



terncurons. and has a recurrent excitation. 2(b) : Related block diagram. 



9 



Parameter Interpretation Value 



A Average excitatory synaptic gain 3.25mV 

B Average inhibitory synaptic gain, slow dendritic inhibition loop 22mV 

C Average inhibitory synaptic gain, fast somatic inhibition loop 20mV 

- Time constant of average excitatory postsynaptic potentials 10ms 
I Time constant of average inhibitory postsynaptic potentials 35ms 

- Time constant of the filter time delay 5ms 
as, a$ Average probability of synaptic contacts in the fast feedback inhibitory loop 0.1 



a? Average probability of synaptic contacts between slow and fast inhibitory neuron 0.8 

Table 2: Parameters interpretations and values of the extended model proposed 
by Wendling and Chauvel (see [50j). The parameters ot\, . . . ,014, J, vq, t and 
Vmax have the same interpretations and values as in Jansen's original model, 
see table Q] 

In this block diagram representation, the PSP functions are given by: 

' h e (t) = Aate~ at 
< hi(t) = Bbte- bt 
h f (t) =Ccte~ ct 

The equations of the extended model hence read: 

yo = 2/5 

y 5 = AaS(y! - y 2 - y 3 ) ~ 2ay 5 ~ a 2 y 

yi = 2/6 

y e = Aa{p(t) + J 2 S(J 1 y )} - 2ay 6 - a 2 y 1 

< 2/2 = V1 (4) 
y 7 =BbJ i S(J 3 y )~2by r -b 2 y 2 

y'3 = 2/8 

y 8 = C cJ 7 S(J 5 y - Jeyi) - 2cy s - c 2 y 3 

2/4 = 2/9 

j/9 = B bS(J 3 y ) - 2by 9 - b 2 y 4 

Numerical values of the parameters This model has been fitted using 
SEEG data, and the authors obtained the values given in the table [2] 

Reduced Wendling-Chauvel's model First of all, one of the most straight- 
forward reduction of the model consists in removing the variables y± and j/9 since 
they are deduced of y 2 and j/7 by the simple formulas: y 2 = J^y^ and 2/7 = J42/9. 
To reduce further the model we proceed in the same way as in Jansen and Rit's 
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case. We make the system dimensionless by introducing the new time t = at 
and proceeding to the change of variables: 



(Y l (T)=Jry l ( T /a) i £ {0, 4} 

Y l {T)=ry l { T /a) i G {1, 2, 3} 

Yi(r) = £ ytir/a) ^{5,9} 

[Yi(r) = 5w(r/o) *G{6, 7, 8} 

We denote by X the interesting signal related to the EEG signal: X — 
Y\ —Y% — Y 3 feeding the pyramidal population and by Z = a§Yo — agY^ = 
a§Yo — a^/aiYa the action of the dendritic inhibitory interncuron population 
on the somatic inhibitory interneuron population. We consider now the variables 
(Y , X, Z, Y 3 , Y 5 ,Y 6 ,Y 7 , Y s ). 

These new variables satisfies the following set of differential equations, where 
we denote for the sake of compactness of notations by a dot the derivative with 
respect to r: 



(5) 



=Y 5 

X =Y 6 ~Y 7 -Y 8 

Y 2 =Y 7 

Y a =Y 8 

Z = a 5 Y 5 - a 6 Y 9 

Y 5 =jS(X)-2Y 5 -Y 

% = ja 2 S( ai Y ) - 2Y 6 -(X + Y 2 + Y 3 ) + P(t) 

Y 7 = jdidonSiaaYo) - 2d x Y 7 - d\Y 2 

% = jd 2 G 2 a 7 S(Z) - 2d 2 Y s - d 2 2 Y 3 

Y 9 =jd 1 G 1 S(a 3 Y )-2d 1 Y 9 -d 2 1 !^^z 

where j, P and S(x) are the same as the one introduced for Jansen and Rit's 
model. We used the notations: 



Gi - -j 



G 2 — § 



(6) 



di = i d 2 = 1 
Using the numerical values in table [2] we obtain 

Gi = 6.76923 G 2 = 6.15385 
di = 0.2857 d 2 = 2 

In this article we are interested in the influence of the other parameters of 
the model together with the input on the fixed points. More precisely, we are 
interested in the codimension two bifurcations of this model with respect to 3 
pairs of parameters: 

1. The input P and the coupling strength j. 

2. The input P and the delays ratio d. 
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3. The input P and the PSP amplitude ratio G. 



3 Influence of the total connectivity parameter 
in Jansen and Rit's model 

We first study the dynamical properties of Jansen and Rit's model. We recall 
in the first subsection the main features described by Grimbcrt and Faugcras in 
[19] . and extend their study to codimcnsion two and three bifurcations. 

3.1 Fixed points and stability 

An interesting property of the system ^ is that the equilibria can be parametrized 
as a function of the state variable X = Y± — Y 2 : 

Y = jS(X) 

Y 2 =a 4 ^jS(a 3 jS(X)) 
a 

P = X + jS(X) - a 2] S{a l3 S{X)) 
The Jacobian matrix at the fixed point is also parametrized by X and reads: 



J(X) = 



( 








1 






















1 


-1 



















1 




-1 


JS'(X) 





-2 










a 1 a2jS'(ja\S(X)) 


-1 


-1 





-2 







\ a 3 a 4 jdS'(ja 3 S(X)) 





-d 2 








-2d 


) 



Although all the dynamics can be parametrized with the variable X, because 
of the complexity of the sigmoidal function, the analytical bifurcation study is 
untractable, and one has to make use of a numerical software in order to solve 
the problerdij. 

In the present case, almost all the calculations can be performed analyti- 
cally in function of the variable X. For this reason, for computing accurately 
the bifurcations of equilibria, we developed our own software implemented using 
Maple® in order to identify our codimcnsion two bifurcations. For simple mod- 
els, the programs we developed give closed-form expressions for the bifurcations 
points. For the continuation of periodic orbits, the analytical study is no more 
possible, and we used both the Matlab® toolbox MatCont [6] and XPPAut 
continuation softwares to study the global bifurcations. The algorithms we use 
and the results we obtain are presented in appendix [A] 

Grimbert and Faugcras used the XPPAut software [9] 
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Figure 3: Projection of a SNIC bifurcation in a subplane of the phase space. 
Three invariant circles are represented: Left. An heteroclinic orbit. Center. 
The two fixed points belonging to the invariant circle merge into a saddle-node 
fixed point. The resulting homoclinic orbit has an infinite period. Right. The 
invariant circle turns into a periodic orbit as the fixed points disappear. 

3.2 Codimension 1 bifurcations 

The dependency of the dynamics to the input firing rate has already been studied 
Grimbert and Faugeras in [19] when parameters are set as in table [T] They 
describe a very rich bifurcation diagram, with the coexistence of two limit cycle, 
one coming from a Hopf bifurcation, and the other collapsing on the fixed points 
manifold, as we show in figure ??. 

The system features two saddle- node and three Hopf bifurcations. One of 
the Hopf bifurcations is subcritical, the other two supercritical. The branch of 
unstable limit cycles originating from this point undergoes a fold bifurcation 
and connects to a family of stable limit cycles of large amplitude that eventu- 
ally collide with a saddle-node bifurcation point and disappear via saddle-node 
homoclinic bifurcation. The periods of these cycles correspond to frequencies 
in the dimensioned model ranging from to 5Hz which is consistent with the 
frequencies of recorded epileptic spikes. It corresponds to what was interpreted 
as epileptic oscillatory activity. 

The other two Hopf bifurcation share the same family of periodic orbits. 
The period of these cycles is almost constant, ranging from 9 to 9.6 in the 
dimcnsionless model, which corresponds in the original model to frequencies in 
the alpha band. 

The system undergoes a saddle-node homoclinic bifurcation (SNIC) at the 
saddle-node bifurcation point, corresponding to a transition between a peri- 
odic orbit and an heteroclinic orbit (see figure [3|), and after this bifurcation, 
the system presents a family of heteroclinic orbits, which are not linked with 
oscillations, since the cycle contains a stable fixed point. 

It had been suggested by Grimbert and Faugeras [19] that this picture was 
quite sensitive to changes in the parameters. They observed that varying any 
parameter by more than 5% resulted in drastic changes in the bifurcation dia- 
gram and the behaviors. This is why we were interested in understanding better 
the appearance of these features and its sensitivity in function of different pa- 
rameters. 
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In this section, we will be particularly interested in the influence of the 
coupling strength j. 

3.3 Effect of the coupling strength and the input current 

Let us first study the bifurcations of the system with respect to the pair of 
parameters (j,P). We first study the bifurcations of equilibria of the system, 
before studying global bifurcations of cycles and the resulting rhythms gener- 
ated. 

We numerically observe that the system undergoes the following bifurcations 
of equilibria (See appendix IA1 and figure ??): 

1. A saddle- node bifurcation manifold (a curve represented in black in figures 
77, Sand©, 

2. An Andronov-Hopf bifurcation manifold (a curve represented in blue for 
subcritical and in red for supercritical bifurcations in the same figures), 

3. A Cusp bifurcation C, 

4. A Bogdanov-Takcns bifurcation BT, 

5. A Bautin bifurcation GH (Generalized Hopf). 

The periodic orbits are also explored using a continuation algorithm. We 
are especially interested in stable cycles which correspond to observable activity 
in the presence of noise. We chose MatCont continuation package developed by 
Kuznetsov, Govaerts and colleagues [SJ [5] and which is very efficient for identi- 
fying bifurcations of periodic orbits. The information of the local bifurcations 
described and the numerical exploration of the bifurcations of limit cycles lead 
us to identify the following generic bifurcations: 

1 . A saddle homoclinic bifurcations linked with the existence of the Bogdanov- 
Takens bifurcation point (plain green curve of figure [S]) . This curve can 
be locally computed using the normal form of the system at this point, 
and continued using a continuation algorithm. Along this curve, a branch 
of limit cycles collapses with the saddle fixed points curve, and the period 
of the related cycles tends to infinity when approaching this curve. When 
numerically continuing this curve we observe that that at the value of 
parameter j related to the Bautin bifurcation, this saddle-homoclinic bi- 
furcation curve collapses with the saddle-node bifurcation manifold giving 
birth to: 

2. A saddle-node homoclinic bifurcation curve ( dashed green curve in Fig. 
[5]) , and connects to a saddle-node limit cycle at a point we note S in 
figure [5] 
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Figure 4: Bautin bifurcation: in region 1 the system has a single stable equi- 
librium an no cycle. Crossing the Hopf bifurcation boundary H_ to region 2 
implies the appearance of a stable limit cycle which survives when we enter 
region 3. Crossing the Hopf boundary H+ creates an extra unstable limit cy- 
cle inside the first one, while the equilibrium regains its stability. Two cycles 
of opposite stability exist in region 3 and collapse at the curve T through a 
fold bifurcation of limit cycles that leaves a single equilibrium. The curves are 
computed in the case of Jansen's model. 

3. A Fold of Limit Cycles exists around the Bautin bifurcation point: a family 
of stable limit cycles and a family of unstable limit cycles collapse and 
disappear along a nondegencrate fold bifurcation of cycles (orange curve 
of [4J. From the Bautin bifurcation, the manifold of FIC is continued and 
we numerically observe that cycles undergo a cusp bifurcation (See figure 
??). The upper branch of limit cycles corresponds to the branch of folds 
of limit cycles generated at the Bautin bifurcation the cycles shrink to a 
single point GH. The lower branch connects to the homoclinic saddle-node 
manifold. At this point, the cycle corresponding to the FIC bifurcation is 
a saddle-node homoclinic cycle (point S). Note that the curve of folds of 
limit cycles originating from the Bautin point can be seen as a function of 
j{P) which is first decreasing then increasing. The point where it changes 
monotony is noted E in figure 0b. 

The full bifurcation diagram is provided in figure [5l The analysis of the 
bifurcation diagram leads to classify the system into 8 classes depending on the 
coupling strength. In each class the system has the same dynamical features 
and the same qualitative behaviors. 
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Figure 5: Codimcnsion two bifurcations of Jansen's model with respect to the 
connection strength j and the input value P. The behavior of the system for 
fixed j can be split into eight zones (A) . . . (H) described in the text. The black 
curve corresponds to the saddle-node bifurcations manifold, the point C the 
cusp bifurcation point, the blue curve corresponds to the subcritical Hopf bifur- 
cations. It is connected to the saddle-node manifold via subcritical Boganov- 
Takcns bifurcation at the point BT. At this point, the saddle-homoclinic bi- 
furcations curve is plotted in green. It exists while j is inferior to the value 
related to the Bautin bifurcation. At this point, the saddle-homoclinic bifur- 
cations curve connects to the saddle-node manifold and generates a curve of 
saddle-node homoclinic bifurcations (dashed green line). The subcritical Hopf 
bifurcation manifold (b;ie curve) is connected to the supercritical Hopf bifurca- 
tion manifold (red curve) through a Bautin bifurcation (point GH). From this 
bifurcation point there is a manifold of folds of limit cycles represented in or- 
ange. The cycles undergo a cusp bifurcation at the CLC point and a saddle-node 
homoclinic bifurcation at the point S. 
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12.38 
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GH 
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3.75 



Table 3: (j, P) coordinates of the different bifurcation and special points of the 
Jansen and Rit's system 
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(a) Equilibria (b) Behaviors 



Figure 6: Case (A): j = 4. (a) equilibria: coordinate Y\ as a function of the 
input current P. For each value of P there exists a unique stable equilibrium, 
(b) behavior of the system. Variation of the X coordinate for different values 
of P: the system always converges to the unique equilibrium. 

3.3.1 Neuro-Computational features 

The different classes of parameters rcsprescnted in color in figure [5] correspond 
to different responses to varying inputs. Eight classes can be distinguishes. We 
observe that the original Jansen and Rit's model is in the zone labelled (D) 
which is of very small extension. 

Non-oscillating behaviors For j < ju2 (zones (A), (B) and (C)), the system 
does not feature any stable oscillations, and therefore the cortical column will 
not oscillate. 

A. For j < j c the system has a unique fixed point and no cycle for any input 
firing rate P. Therefore, when the input firing rate is fixed, for any initial 
condition, the activity of the column will converge towards this unique 
equilibrium. In that case the cortical column has a quite trivial behavior: 
it has no oscillatory activity and converges to rest whatever the input. 

B. For jc < j < 3BT-, the system undergoes two saddle- node bifurcations 
when varying the parameter P. Depending on the value of P, the sys- 
tem has one, two or three fixed points and no cycle (see figure [7]). For 
P $ [P\,P2\ (i.e. not between the values of the the two saddle-node bi- 
furcations), there is a unique fixed point which is stable and the system 
converges to this fixed point. When P G [P-^^Pi), there are three fixed 
point, one is unstable and the other two stable. The system is therefore 
bistable: depending on the initial condition, the activity will converge to- 
wards one or the other stable fixed point, corresponding to an up-state 
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Figure 7: Case (B): j = 8. (a) equilibria (coordinate Y\ as a function of the 
input current P), and stability. Continuous line: stable equilibrium, dashed 
line: unstable. For P [Pi,p2] there exists a unique stable equilibrium, and 
for P G (Pi,P 2 ) there exist three equilibria, two stable and one unstable, (b): 
Behavior of the system for different input and initial conditions. Purple: P = 
— 10 and green: P = 10 : the system always converges to the unique equilibrium. 
Orange: P = 2: bistability. We observe that the activity resembles a real evoked 
potential. 



and down-state activity. The system also presents hysteresis when contin- 
uously varying the input in this zone of inputs. Eventually, it can switch 
between the two stable fixed point if perturbed. 

C. For jsT < 3 < 3hi- The system has two saddle-node bifurcations and a 
subcritical Hopf bifurcation (see figure |5J). Therefore, the system has a 
unique stable fixed point when the input P is not between the two saddle- 
node bifurcation points (i.e. P [Pi,iV]), and the system converges 
towards this fixed point. For P between the first saddle-node bifurcation 
value and the Hopf bifurcation (i.e. P € (Pi,Ph)), the system has 2 
unstable fixed points and a stable fixed point, and generically converges 
towards the stable fixed point. In the zone between the Hopf bifurcation 
and saddle homoclinic point (P € (Pff f Psfe)) , the system presents two 
stable fixed points, an unstable fixed point and an unstable limit cycle. 
Depending on the initial condition, the system will either converge to one 
or the other fixed point. When P is greater than the saddle- homoclinic bi- 
furcation value and below the greatest saddle-node value (P G (Psh, P2)), 
the system has two stable fixed points and an unstable fixed point, and 
its behavior is similar to the behavior in the previous case. We can see 
that the system returns to equilibrium via oscillations. Hence in the pres- 
ence of noise, the system will present oscillations at a certain frequency 
superimposed to its noisy behavior. 
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(a) Equilibria 



(b) Behaviors 



Figure 8: Case (C): j = 11. (a): Equilibria and bifurcations (coordinate Y\ as a 
function of the input current P), and stability. Blue lines: Continuous : stable 
equilibrium, dashed : unstable. Dashed purple line: unstable cycles. Green 
cycle: saddle homoclinic orbit. For P [P\,P2\ there exists a unique stable 
equilibrium, and for P G {P\ 1 P2) there exist three equilibria: for P G (Pi,Ph) 
two unstable and one stable, and no limit cycle. For P G (Ph,P2), two stable 
and one unstable. For P G [Ph,Psii\, there exists an unstable limit cycle, 
(b): Behavior of the system for different input and initial conditions. Purple: 
P = — 10 and green: P = 10 : the system always converges to the unique 
equilibrium. Orange: P = 0: bistability and damped oscillations. We observe 
that the activity resembles a real evoked potential. 
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Rhythmic activity For values of the connectivity greater than jjj 2 , the sys- 
tem will always present a supercritical Hopf bifurcation and therefore a stable 
periodic orbit, corresponding to a rhythmic activity of the column. 



D. For j € [7j? 2 ,Jb], the system undergoes two saddlc-nodc bifurcations, two 
supercritical and a subcritical Hopf bifurcations, and one fold of limit 
cycles. It is the case of the original Jansen and Rit's model [TJ We label 
the bifurcation points as in figure [51 

For P < Pi , the system has a unique stable fixed point and for any initial 
condition, it converges towards this point. 

For P\ < P < Ph 1 , the system has three fixed points, one of which is 
stable and the other two unstable, and the system still converges towards 
the unique fixed point. 

For Pjj 1 < P < Ph 2 -, the largest fixed point becomes stable and there 
is an unstable limit cycle (represented in dashed purple in figure [9]) . In 
this region, the system converges towards one of the two stable fixed point 
depending on the initial condition of the system. 

For Pjj 2 < P < P2, the largest fixed point loses stability via a subcriti- 
cal Hopf bifurcation and a stable periodic exists with frequency of about 
10Hz corresponding to purely alpha activity. In this region, the system 
either oscillates around the periodic orbit and presents alpha-activity, or 
converges to the stable fixed point. 

For P 2 < P < Pflc, the system has no fixed point and two stable limit 
cycles: the cycle corresponding to the continuation of the subcritical Hopf 
bifurcation which corresponds to alpha activity, and a large amplitude 
cycle with a frequency ranging from to 5 Hz which corresponds to an 
epileptiform activity. The system selects one of these cycles depending on 
the initial condition, and can switch from one activity to another when 
the system is perturbed. Assume that we slowly increase the input P. 
If the system was in the down equilibrium state, it will converge to the 
epileptic cycle when crossing the bifurcation, and if it was in the up state 
to the alpha cycle. 

For Pflc < P < Ph 3 , the system has a unique stable trajectory which is 
a periodic orbit with frequency close to 10Hz, and for any initial condition 
the system will converge towards this cycle. Eventually, for P > Ph 3 the 
system has a unique fixed point and for any initial condition, the state 
converges towards this fixed point. 

E. For J G [jEijcii], the number and stability of the fixed points are the 
same as in the previous case. But in this case, the structure of the cycles 
is more complex. The family of unstable periodic orbits originating from 
the subcritical Hopf bifurcation is connected to the family of limit cycles 
of the supercritical Hopf bifurcation H2 associated with the smaller P 
value, and the branch of limit cycles of the supercritical Hopf bifurcation 
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(a) Equilibria 



(b) Behaviors 



(c) Periods 



Figure 9: Case (D): j = 12.285. (a): Equilibria and bifurcations (coordinate 
Y\ as a function of the input current P), and stability. Blue lines: Continuous: 
stable equilibrium, dashed: unstable. Purple lines: dashed: unstable cycles, 
continuous: stable. Green cycle: saddle homoclinic orbit, orange cycle: fold of 
limit cycles (see text) (b): Behavior of the system for different input and initial 
conditions. Purple: P = — 10 and green: P = 10: the system always converges 
to the unique equilibrium. Red and orange: P = 2.3: bistability of cycles: 
orange: epileptic spikes, red: alpha activity, (c): period of the limit cycles. 
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H 3 associated with the greatest P value undergoes two fold bifurcations 
of cycles and disappears via saddle-node homoclinic bifurcation (see fig- 
ure [TCI and labels herein). For P < Pflc\ the behavior of the system is 
exactly the same as in case (D): the system generically converges to the 
unique stable fixed point for P < Ph 2 , and either converges to the stable 
fixed point or to the stable cycle depending on the initial condition. For 
PflCx < P < Pflc 2 , the system return to the down-state equilibrium. 
For Pflc 2 < P < P2, the system has a stable fixed point (down-state 
equilibrium), a stable limit cycle corresponding to alpha-like activity and 
an unstable limit cycle. In this zone of input the system will experience 
bistability between alpha-activity and rest. For P 2 < P < Pflc 3 , the 
system has no stable fixed points and three cycles, one of which is un- 
stable, another one corresponding to alpha activity and the third one to 
epileptiform activity. For P > Pflc 3 the behavior is the same as in the 
case (D) for P > Pflc- while P < Ph 3 the system presents purely alpha 
oscillations, and for P > Ph 3 the system converges to the unique up-state 
equilibrium. 

F. For jai < j < 3Hn the system has two subcritical Hopf bifurcations 
whose family of limit cycles are connected, and a supercritical Hopf bi- 
furcation whose limit cycles undergo two saddle-node bifurcations of limit 
cycles and collapse on the saddle-node manifold via saddle-node homo- 
clinic bifurcation (see figure ITTj) . For P < Ph 1 the system converges to 
the down-state equilibrium. For Ph 1 < P < Ph 2 , the system is bistable 
and either converges to the upstate equilibrium or to the downstate equi- 
librium depending on the initial condition. For Ph 2 < P < Pi the system 
converges to the downstate equilibrium. Therefore for P < Ph 2 , the sys- 
tem only presents damped subthreshold oscillation and no real rhythmic 
activity. For Pi < P < PflCi the system is in a pure epileptic activity, 
for PflCx < P < Pflc 2 the system presents both epileptic activity and 
alpha activity, and for Pflc 2 < P < Ph 3 the system only presents alpha 
activity. In this case, when slowly varying the input P, the system will 
always present epileptic activity: indeed, it is the only stable activity in a 
certain range of parameters. 

G. For Jh 1 < j < jcLC, the system has two saddle-node bifurcations and 
a supercritical Hopf bifurcation, whose limit cycles undergo two saddle- 
node bifurcations and disappear via a saddle-node homoclinic bifurcation. 
For all P < P 2 , the system always converges to the downstate equilib- 
rium. For P 2 < P < PflCi the system presents epileptic spikes, and for 
PflCi < P < Pflc 2 bistability with epileptic spikes and alpha activity. 
For Pflc 2 < P < Phx the system presents only a stable alpha activity 
and for P > Ph 1 the system always returns to an upstate equilibrium 
whatever the initial condition (see figure [T2|) 

H. For j > jcLC the system features two saddle-node bifurcations and a 
supercritical Hopf bifurcation whose family of limit cycle is regular and 
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(a) Equilibria 



(b) Periods 



Figure 10: Case (E): j = 12.42. (a): Equilibria and bifurcations (coordinate 
Y\ as a function of the input current P), and stability. Blue lines: Continuous: 
stable equilibrium, dashed: unstable. Purple line: dashed: unstable cycles, 
plain: stable. Green cycle: saddle homoclinic orbit, orange cycles: folds of limit 
cycles (see text), (b): Behavior of the system for different inputs and initial 
conditions. Purple: P = —4 and green: P = 10: the system always converges 
to the unique equilibrium. Red and orange: P = 2.4: alpha and epileptic 
activity, light blue: P = 1.1: rest and alpha oscillations, pink: P = 4: purely 
alpha activity. E and time removed because arxiv refuses large files. 
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(a) Equilibria 




(b) Behaviors 



(c) Periods 



Figure 11: Case (F): j = 12.5. (a): Equilibria and bifurcations (coordinate Y± 
as a function of the input current P), and stability. Blue lines: Continuous: 
stable equilibrium, dashed: unstable. Purple line: dashed: unstable cycles, 
continuous: stable. Green cycle: saddle homoclinic orbit, orange cycles: folds 
of limit cycles (see text), (b): Behavior of the system for different inputs and 
initial conditions. Purple: P = —2 and green: P = 9 : the system always 
converges to the unique equilibrium. Dark blue: P = 2.05: purely epileptic 
activity: slow waves of high amplitude. Red and orange: P = 2.3: alpha and 
epileptic activity coexist, light blue: i"24= 5: only alpha oscillations. 




Figure 12: Case (G): j = 12.7. (a): Equilibria and bifurcations (coordinate Y\ as 
a function of the input current P), and stability. Blue lines: Continuous: stable 
equilibrium, dashed: unstable. Purple line: dashed: unstable cycles, continuous: 
stable. Green cycle: saddle homoclinic orbit, orange cycles: folds of limit cycles 
(see text), (b): Behavior of the system for different inputs and initial conditions. 
Purple: P = and green: P = 11: the system always converges to the unique 
equilibrium. Orange: P — 2.3: purely epileptic activity: slow waves of high 
amplitude. Pink: P = 3.1: alpha and epileptic activities coexist, light blue: 
P = 5: only alpha oscillations. 
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disappears via saddle- node homoclinic bifurcation (see figure [12]). For 
P < P2 or P > Phj the system converges to the unique fixed point, and 
for P2 < P < Ph! the system presents oscillations. This case is very 
interesting from a neuro-computational point of view. Indeed, the family 
of limit cycles created presents a stable region of oscillations around 10Hz 
corresponding to alpha activity as in the previous case, for a quite large 
set of inputs P a < P < Ph x ■ At P = P a it abruptly switches from alpha 
activity to theta activity where it stays for Pg < P < P a and eventually 
switches regularly to delta activity for P 2 < P < P$. In this region of 
parameter therefore the system presents the rhythm of the normal sleep 
activity. 
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Figure 13: Case (H): j = 14. (a) Equilibria and bifurcations (coordinate Y\ 
as a function of the input current P), and stability. Blue lines: Continuous: 
stable equilibrium, dashed: unstable. Purple line: dashed: unstable cycles, 
plain: stable (see text), (b) Behavior of the system for different inputs and 
initial conditions. Up: P = 10: alpha activity. Middle: P = 5 theta activity. 
Down P = 2.3: delta activity. 27 



4 Influence of other parameters in Jansen and 
Rit's model 



4.1 Effect of the PSP amplitude ratio G 

The bifurcation structure as a function of the PSP amplitude ratio G is very 
similar to the one corresponding to the coupling strength j. The system also 
presents a cusp, a Bogdanov-Takens and a Bautin bifurcation, two branches of 
saddle-node of limit cycles collapsing at a cusp of limit cycles. The same types 
of behaviors are observed. We decompose here again the bifurcation diagram 
into zones depending of G, and use the same notations as in the previous case. 
The bifurcation diagram and the decomposition in zones is given in figure [Til 
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Table 4: (G, P) coordinates of the different bifurcation and special points of the 
Jansen and Rit's system 
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Figure 14: Codimention 2 bifurcations in Jansen and Rit's model with respect 
to the PSP amplitude ratio G and the input P. The behavior of the system 
for fixed values of G can be split into eight zones (A) . . . (H) described in the 
text. The black curve corresponds to the saddle-node bifurcations manifold, the 
point C to the cusp bifurcation point, the blue curve to the subcritical Hopf 
bifurcations manifold. It is connected to the saddle-node manifold via a sub- 
critical Boganov-Takens bifurcation at the point BT. At this point, the saddle- 
homoclinic bifurcations curve is plotted in green and connects to the saddle-node 
manifold via saddle-homoclinic bifurcation, becoming a saddle-node homoclinic 
bifurcation curve (dashed green line) . The subcritical Hopf bifurcation manifold 
is connected to the supercritical Hopf bifurcation manifold (red curve) through 
a Bautin bifurcation (point GH). From this bifurcation point there is a manifold 
of folds of limit cycles represented in orange. Cycles undergo a cusp bifurcation 
at the CLC point and a saddle-node homoclinic bifurcation at the point S. (See 
values of the parameters in table 2]) 

The behaviors in the eight zones are very similar. The behavior in the zone 
(A) is similar to the behavior of the system in the zone (B) in figure [3 zone (B) 
to zone (C) of figurc[SJ zone (C) to zone (D) of figure [SI zone (D) to (E), (E) to 
(F) and (F) to (G) and (G) to (H). The only new kind of behavior corresponds 
to zone (H) where there exist no saddle node bifurcation but a supercritical 
Hopf bifurcation which corresponds to the existence of a cycle. 
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Figure 15: Codimcnsion 2 bifurcations in Jansen and Rit's model with respect 
to the delay ratio d and the input parameter P. 

4.2 Effect of the delay ratio d 

The delay ratio d features the same bifurcation structure as in the previous 
two cases with a cusp, a Bogdanov-Takens, a Bautin and a cusp of limit cy- 
cles bifurcations, together with saddle-homoclinic and saddle-node homoclinic 
bifurcations. The bifurcation diagram is shown in figure [TBI 

4.3 Sensitivity to the connection probability parameters 

The structure of cycles we studied is unaffected by changes in a\, 03 and ct\. 
For the parameter 012, the picture is slightly different. Indeed, the bifurcation 
diagram presents a codimension three bifurcation corresponding to the cusp case 
of the degenerate Bogdanov-Takens bifurcation. With respect to these parame- 
ters, the cusp bifurcation point is also a point of a Andronov-Hopf bifurcation. 
Note that at this very point, the bifurcation diagram is very close to the one of 
Wcndling and Chauvel's model (see section[5]) and therefore will locally generate 
the same behaviors as those of this more complex model. 

This codimcnsion three bifurcation corresponds to the values DBT : {P = 
3.236, a 2 = 0.365}. 
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(a) Codim. 2 bifurcations (b) Codim. 2 bifurcations 
with respect to 03 and P with respect to 04 and P 
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(c) Codim. 2 bifurcations with respect to 02 and P 



Figure 16: Codimcnsion 2 bifurcations of Jansen and Rit's model with respect to 
the connection probability parameters a% (removed because Arxiv refuses large 
files)., i = 1, • ■ • ,4. In the first three cases no new bifurcations appear: there are 
only a saddle-node and a Hopf bifurcation manifolds, with possibly Bogdanov- 
Takens and Bautin bifurcation points. In case (d) we observe a codimcnsion 
three degenerate Bogdanov Takens (DBT) bifurcation point. 
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5 Analysis of Wendling and Chauvel's model 



We now turn to the study of the Wendling and Chauvel's model of hippocampus 
activity [50j . As stated, this model was initially introduced in order to be 
able to reproduce the kind of fast activity encountered at the onset of seizures 
in limbic structures. From many experimental observations, the authors take 
into account slow and fast currents, which resulted in the model presented in 
section 12.21 We show here that codimension one and two bifurcations with 
the standard parameters for biologically plausible processes do not lead to the 
emergence of this fast activity. In order to reproduce this kind of activity 
Wendling and Chauvel used changes in the amplitude of the PSPs. We present 
the mathematical emergence of this fast activity via a bifurcation analysis, and 
discuss the effect of the noise. 

5.1 Fixed points of the model 

Wendling and Chauvel's dynamical system is even more intricate than Jansen 
and Rit's, since it has ten dimensions, sigmoidal nonlincaritics and component 
mixing. Nevertheless in this case again, the fixed points of the system can be 
parametrized as a function of the state variable X — Y\ — Y 2 . This manifold is 
given by the following set equations: 

'Y =jS(X) 
Y 2 = jG^± S{adS{X )) 

Z = -2G^ S (a 3 jS(X)) + a 5 jS(X) 

< Y 3 = 2G^L S (-^S(a 3 jS(X)) + a 5 jS(X)) 

P =X + Y 2 (X)+Y 3 {X)-a 2 jS(a 1J S(X)) 

Yi =0 Vie {5... 9} 

The input firing rate P is an important parameter. This is why we always 
consider bifurcations with respect to this parameter. 

5.2 Bifurcations of the Wendling and Chauvel model 

5.2.1 A degenerate codimension three Bogdanov-Takens bifurcation 

We first consider, similarly to Jansen and Rit's case, the codimension two bi- 
furcations with respect to the input P and the total connectivity parameter 
3- 

We numerically observe that for the original values of Wendling and Chau- 
vel's model, the system features a saddle-node and a Hopf bifurcation manifolds, 
a cusp, Bautin a codimension three bifurcation: a degenerate Bogdanov-Takens 
bifurcation. Figure [17] shows the corresponding bifurcation diagram which is 
split into four zones labelled from A to D. 
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Figure 17: Bifurcations of the Wendling and Chauvel's model with respect to 
the parameters P and j. The other parameters are given by the values in table 
[2] We observe a saddle-node bifurcation manifold (plain black curve), a Hopf 
bifurcation manifold (plain blue curve), a fold bifurcation limit cycles manifold 
(plain orange curve), a Bautin bifurcation, and a codimcnsion three degenerate 
Bogdanov-Takens bifurcation (DBT). 
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Table 5: (j, P) coordinates of the different bifurcation and special points of the 
Wendling and Chauvel's model. 
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Figure 18: Wcndling-Chauvel's model in case (A): unconditional convergence 
to a unique fixed point, j=5, blue P=-5, green: P=10. 

A. In this zone, the system does not present any bifurcation nor any cycle, 
and the system returns to equilibrium for any initial condition and any 
input rate. 

B. In this zone, the system features two saddle-node bifurcations and a sub- 
critical Hopf bifurcation generating unstable limit cycles. The family of 
cycles undergoes a fold of limit cycles, and the system displays slow oscil- 
lations of large amplitude, corresponding to epileptic-like activity. 

C. In this case, the system presents alpha activity together with epileptic 
spikes of low frequency and high amplitude because of the presence of two 
fold bifurcations of limit cycles. 

D. Case (D) corresponds to the case where the system goes continuously from 
alpha activity to slow waves that can be interpreted as sleep waves. 

We observe that, so far, i.e. in the case where the parameters are set using 
biophysical values, Wcndling and Chauvel's system only features behaviors that 
are also featured by Jansen and flit's model. In particular it does not display 
the kind of fast activity that motivated its development. 

5.3 Behaviors of Wendling and Chauvel's model and bi- 
furcations 

In [50], the authors provide a partition of the parameter space (A,B,C) of 
the original system corresponding to different behaviors of the system, and 
distinguish between normal activity, narrowband activity, sporadic spikes and 
fast activity. The partition has been computed simulating different EEG traces 
for a Brownian input rate. The reduction of the model proposed in section [2. 2 1 
shows that in the space (A, B,C), bifurcations can be described via a bifurcation 
analysis in the space (Gi = B/A,G2 = C/A). Nevertheless, one has to be 
careful that the coefficient A influences the input firing rate since P = j^P- 
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(a) Case (B): Bifurcation diagram 




(b) Case (B): Periods of the limit cy- (c) Case (B): Time trajectories 
cles 



Figure 19: Wendlmg-Chauvel's model in case (B): Epileptic Spikes. Parameters: 
j=8. (c): blue: P=2, purple: P=3.5, orange P=4, green: P=5.5, cyan P=7, 
black P=9. 
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(a) Case (C): Bifurcation diagram 




(b) Case (C): Periods of the limit cycles 




(c) Case (C): Time trajectories 



Figure 20: Wendling-Chauvel's model in case (C): alpha and epileptic waves. 
Blue and red: P=9.7, pink P=7, orange P=5. 
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(a) Case (D): Bifurcation diagram (b) Case (D): Periods of the limit cy- 
cles 




(c) Case (D): Time trajectories 



Figure 21: Wendling-ChauvePs model in case (D): Sleep waves, j=13, P=5, 15 
and 22. 
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In the case where A is small, no oscillatory activity is seen in the region of 
interest of the parameters: the system is quite far from any bifurcation. For 
slighter higher values of this parameter, the system presents two Hopf bifur- 
cations which share the same family of periodic orbits leading to narrowband 
activity (a, 8), which corresponds to the case A = 3 in [50], as represented in 



figures 22(a) and 22(b) The fast activity (/3, 7) observed numerically has no 
deterministic equivalent. It corresponds to a purely stochastic phenomenon: 
there exists a unique stable fixed point in the region presented by Wendling and 
Chauvel in [50], and the Jacobian matrix has complex eigenvalues. The noisy 
inputs destabilize the system, which moves randomly around the fixed point. 
The way the system returns to equilibrium corresponds to oscillations at the 
frequency corresponding to the inverse of the imaginary part of the eigenvalue. 
This is why the spectrum of the signal generated in the neighboorhood of this 
fixed point presents a peak in the Fourier spectrum, but this spectrum is quite 
spread out and the amplitude of the signal is quite low, resembling purely noisy 
activity. 

For larger values of A, folds of limit cycle appear. In that case, we observe 
a small region of narrowband activity (a, 0), immediately followed by rhythmic 
epileptic spikes corresponding to a family of limit cycles of large amplitude and 
low frequency. Sporadic spikes appear in the region where an unstable cycle 
exists, which is coherent with the interpretation given in the case of Jansen and 
Rit's model. We also see that it corresponds to a random attraction to a nearby 
epileptic cycle: sporadic spikes in Wendling and Chauvel's model appear only for 
regions close to the epileptic activity. Around the Hopf bifurcation, the system 
features a stable fixed point whose Jacobian matrix has non-real eigenvalues. 
The imaginary part of the eigenvalue is close to the value corresponding to a/ 9 
activity near the Hopf bifurcation, and then increases, which corresponds to a 
transition between narrowband and fast activities. This picture remains valid 
when increasing further the parameter A, as observed in |50j . 

For intermediate values of the parameter A, we observe deterministic fast 
activity on a quite strange manifold, as presented in figure [2"31 existing for quite 
a small range of values of the parameter A. But the phenomenology stays the 
same even when deterministic cycles do not exist, as already explained. Note 
that these cycles do not appear in Wendling and Chauvel's model with the 
original parameters, since these cycles and these types of rapid period do not 
appear in the bifurcation diagram in figure RTFl 
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(a) Small values of A: 



Cycles 



(b) Small values of A: Periods 
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Figure 22: Behaviours and bifurcations in Wendling and Chauvel's model. 
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(a) Codimension two bifurcation diagram 



(b) Period of the cycles 



Figure 23: (a): Periods of the limit cycles. Purple band: Very fast waves, Pink: 
fast activity 7), blue: narrowband activity (a, 6), orange: slow activity, (b) 
Codimension two bifurcation diagram with respect to G\ and G2. We observe 
numerically the emergence of the different behaviors observed by Wendling and 
Chauvel. (c) Shape of the cycles in the space (Gi,Yq, X). 



6 Applications: Epileptic seizures dynamics 

Interictal spikes appear at a saddle-node homoclinic bifurcation: for inputs lower 
than the bifurcation value, the system shows a regular behavior corresponding 
to the existence of a stable fixed point. Perturbing this equilibrium by taking 
into account random activity leads to the appearance of non-rhythmic spikes 
which can be interpreted as interictal spikes. 

In the cases (C) to (E) of the diagram in figure [3 interictal spikes appear, 
but do not always correspond to a pre-ictal activity. Indeed, the destabilization 
of the stable fixed point may lead the system around the epilepsy limit cycle, but 
since this cycle is not the only stable behavior, the apparition of random spikes 
will hardly lead to a seizure. Indeed, it appears that the attraction basin of the 
epileptic cycles is quite small. When destabilizing the fixed point, we are very 
likely to present a spike, because at the saddle-node bifurcation point is also 
a point of the saddle-node homoclinic cycle, part of the branch of limit cycles 
corresponding to epileptic spikes. But the noise will soon make the solution quit 
the attraction basin of the cycle and return to the other stable state. In this 
case, isolated random spikes correspond to some kind of transition from rest to 



oscillatory activity (see figure 24(a) ). In this case we can observe the emergence 
of paroxismal depolarizing shifts, as described for instance in [3"Tl chapter 46]. 

In the cases (F) and (G), epilepsy is the only possible behavior in a certain 
range of input parameters. Therefore, interictal spikes really predict a move 
towards the epileptic seizure, and therefore are an indication of a transition to 
a seizure. The type of noise input chosen, the standard deviation of this noise 
and the equilibrium perturbed impact the structure and the properties of the 
spike train generated, but the qualitative behavior is always the same. When 
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the inputs keep the system reasonably far from the saddle-node homoclinic 
bifurcation, spikes have a very small probability to occur and in practice are 
hardly observed. When getting closer to the bifurcation point, isolated spikes 
appear, sec figure [24(b)] Their occurence depends on the value of the input firing 
rate and the noise standard deviation. When increasing further the mean input 
and getting closer to the bifurcation points, spikes become almost rhythmic 
and several superimposed action potentials can appear, which are known to 
be related with paroxismal depolarizing shifts (PDS), see figure [24(b)] When 
increasing again the stimulation, we are in an epileptic crisis, and we observe 
rhythmic high-amplitude epileptic spikes. 




(a) j=12.285 



(b) j=12.7 



Figure 24: Different behaviours in Jansen's reduced model driven by a Brownian 
noise P of mean /x = 1.8 and standard deviation of .5. (a) Interictal spikes 
leading to oscillatory activity. The parameters are those given in (b) 
Interictal spikes and paroxismal depolarizing shift (PDS). The parameters are 
those given in ^ except for j = 12.7. 



7 A possibly pathological scheme 

It has been observed by many experimentalists that an inherited epilepsy could 
be linked with an hyperinnervation (see e.g. [?, ?, ?]). This would correspond 
in our models to an increase of the value of the total connectivity j (as well as 
possibly some changes in the proportions coefficients 014). As a first approxima- 
tion, we observe that small but pathological increases of the total connectivity 
coefficient j can lead to the "epileptic parameter zones" (F) and (G). Let us 
consider that in the resting state the column receives inputs corresponding to 
a normal (fixed point) behavior, and that some kind of stimulus results in a 
slow increase of the mean level of inputs. In that case, the column first behaves 
normally, then goes through a slow transition to epileptic activity during which 
interictal spikes and paroxismal depolarizing shifts are likely to occur (see figure 

USD- 



41 




200 400 600 SOC 1000 1200 1400 1600 1800 



Figure 25: Seizure in the pathological scheme. We consider an small hypcrinner- 
vation corresponding to j = 12.7. The inputs are considered to be random, with 
a slowly increasing mean \i = 1.5 + 10 — 3 t and constant standard deviation of 0.4. 
The purple zone corresponds to normal activity, the pink zone to the onset of 
the seizure: we observe randomly distributed spikes and paroxismal depolariz- 
ing shifts. The blue zone corresponds to the seizure itself: we observe rhythmic 
activity with high amplitude oscillations. The orange zone corresponds to the 
end of the seizure, with a return to alpha activity. We observe that the signal 
naturally waxes and wanes. 

8 Mathematical interpretation of EEG record- 
ings 

What makes a good neural mass model? This question is very far reaching 
and the answer clearly depends on the phenomena one wants to account for. 
From the mathematical point of view, we have been able in the different cases 
studied in this paper to interpret oscillatory behaviors via a bifurcation analysis 
of cycles. Lennaert van Veen and David Lilcy have studied another model of 
EEG signal in [47] . They observed a bifurcation diagram very similar to the one 
of Janscn and Rit's model but in their case the saddle connection S corresponds 
to a Shilnikov Saddle-Node bifurcation [33] . 

The mathematical approach based on bifurcation theory that we developed 
in order to understand the behavior of neural masses was already successfully 
applied to single neuron modeling. In this domain, bifurcations have been re- 
lated to the excitability properties of the underlying dynamical system since the 
pioneering paper of Rinzel and Ermentrout [37] in the late 80s, and the canonical 
models studied by Ermentrout and Kopell. This approach was then generalized 
to the study of the excitability properties of the neurons (see, e.g., [57] [2S] 
for reviews). This approach naturally led to introduce such phenomcnogical 
models of neurons as the nonlinear twodimcnsional intcgrate-and-fire neuron 
[25] [2] [40l [41] . In the case of neural masses there are two outstanding questions 
that we need to answer: first, can we relate EEG behaviors with bifurcations 
of the underlying dynamical system? and second, can we propose a "minimal 
model" that is able to reproduce phenomenologically EEG signals? We next 
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address these two questions. 

8.1 Bifurcations and oscillatory behaviors 

EEG signals are characterized by different phenomena which can be linked to 
mathematical objects: 

• The normal activity is related to the existence of a stable fixed point. 

• Oscillatory activity, whatever its frequency band, is related to the exis- 
tence of a family of limit cycles. 

• Fast activity can be either related to the existence of a family of limit 
cycles having a high frequency, or to the dcstabilization through noise of 
a stable equilibrium whose Jacobian matrix has complex eigenvalues. 

• Epileptic spikes correspond to low frequency, large amplitude, oscillations. 
In general these oscillations appear suddenly when varying the parameters 
which corresponds to the advent of a seizure. Low frequency oscillations 
can be related to the existence of homoclinic orbits or homoclinic bifurca- 
tions which are for instance caused by the presence of a Bogdanov-Takens 
bifurcation. The simplest mathematical object accounting for a sudden 
appearance of large amplitude oscillations is the fold bifurcation of limit 
cycles. 

• In some cases it can be interesting for the model to feature a bistabil- 
ity between an epileptic activity and a rhythmic activity. The simplest 
bifurcation accounting for this coexistence is the cusp of limit cycles. 

Note that all these mathematical objects are present in the two models that are 
studied in this paper. 

8.2 A minimal model? 

A suitable model for reproducing EEG activity could therefore contain the bir- 
furcations just cited. In order to obtain a simple model of this kind, we build 
on the observation that both Jansen and Rit's and Wendling and Chauvel's 
models feature a codimension three bifurcation with the universal unfolding of 
the degenerate Bogdanov-Takens bifurcation. This bifurcation corresponds to 
the cusp case of the singularity of planar vector fields with nillpotent Jacobian 
matrix in the cusp case. This bifurcation was studied by Dumortier, Roussarie 
and Sotomayor in [7] , and then complemented with the elliptic, saddle and focus 
cases by Dumortier and collaborators in [5]. 

The center manifold of this bifurcation is two-dimensional and the normal 
form of this bifurcation reads on the center manifold: 




(7) 
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This bifurcation is generic, as proved in [5], therefore any system that can 
be put in this form with higher order terms will present the same bifurcation 
diagram. Dumortier and collaborators proved also that locally around the De- 
generate Bogdanov-Takens point (DBT), the system presents a curve of cusp 
bifurcations, Bogdanov-Takens bifurcations, two saddle-node bifurcation curves, 
a Hopf bifurcation curve, a curve of fold bifurcation of limit cycles, a Bautin 
(generalized Hopf) bifurcation curve and a saddle-node homoclinic curve. 

Hence, when varying the input rate and the connectivity strength, this model 
features exactly the same bifurcations as Wendling and Chauvel's with their 
original parameters. 

9 Conclusion 

We have presented in this article some results arising from the study of bifur- 
cations in neural mass models in order to understand the origin of some brain 
rhythms and epileptic seizures. The study of the different cycles of the system is 
of great interest. The effect of the different parameters on the model behaviour 
can be understood in the light of these bifurcations. The effect of noise on such 
systems, and the relations between the dynamics of these deterministic models 
and the related mean field equation is a noble endeavor which will probably 
provide great insights on the effect of noise in the apparition of seizures and in 
the rhythms of the brain. 

A A symbolic and numerical bifurcation algo- 
rithm 

In this appendix we present an algorithm to compute formally or numerically 
local bifurcations for vector fields implemented in Maple®. This algorithm is 
based on the closed form formulations for genericity and transversality condi- 
tions given in textbooks such as [201 133] - For our numerical experiments, we 
implemented a precise and efficient solver of equations based on dichotomy. 
This algorithm controls the precision of the solution we are searching for, and is 
way faster than the native f solve Maple application. It is the one we used in 
this article to study the Jansen and Rit's and Wendling and Chauvel's models 
bifurcations. 

A.l Saddle-node bifurcation manifold 

We recall that given a dynamical system of the type x = f(n,x) for x £ R" 
and € R, the systems undergoes a saddlc-nodc bifurcation at the equilibrium 
x = Xo , /i = fiQ if and only if (see e.g. [20l Theorem 3.4.1.]): 

(SN1). D x f{fj,Q,Xo) has a simple eigenvalue. Denote by v (resp w) the right 
(resp. left) eigenvector. 
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(SN2). (w,df/dfi)(x ,(x ))^0 

(SN3). {w,(Dlf{fio,x ))(v,v)}^0 

The algorithm we use to identify the saddle-node bifurcation manifold is a 
straightforward application of this theorem, and consists in: 

• solving the implicit equation det(D x f(po,xo)) = on the fixed-points 
manifold. We obtain the equilibria where the Jacobian determinant van- 
ishes, i.e. where there is a null eigenvalue. 

• We then compute the left and right eigenvectors of the Jacobian matrix 
at these points: 

— check that the dimension of the eigenspace related to the eigenvalue 
is 1 



- check conditions (SN2) and (SN3) 



At the points where the determinant of the Jacobian matrix vanishes and 
the conditions are not satisfied, we are led to consider the two classical cases of 
the cusp and the Bogdanov-Takens bifurcations. 



A. 2 Cusp bifurcation 

At a cusp bifurcation point, the Jacobian matrix of the system has a null eigen- 
value and the first coefficient of the normal form (given by condition (SN3) ) 
vanishes. In that case, under the differential transvcrsality and gencricity con- 
ditions of [331 lemma 8.1], a smooth change of coordinates will put locally the 
system in the normal form of the cusp bifurcation. Our algorithm numerically 
checks the two conditions at these singular points. 



A. 3 Bogdanov-Takens bifurcation 

If along the saddle node manifold a second eigenvalue vanishes, under the gener- 
icity and transversality conditions of chapter 7.3], a smooth change of coor- 
dinates will locally put the system in the normal form of the Bogdanov-Takens 
bifurcation. These conditions are also numerically checked by our algorithm. 



A. 4 Andronov-Hopf bifurcation manifold 

Changes in the stability of fixed points can also occur via Andronov-Hopf bi- 
furcations. In this case, the real part of an eigenvalue crosses but not its 
imaginary part. Theoretically, the dynamical system x = f(x,n) undergoes a 
Hopf bifurcation at the point x = xq, \i = fiQ if and only if (see [20\ Theorem 
3.4.2]): 

(HI). D x f(xo, no) has a simple pair of pure imaginary eigenvalues and no other 
eigenvalues with zero real part. Denote by A(/z) the eigenvalue which is 
purely imaginary at fio- 
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(H2). Zi(xq,/xo) 7^ where l\ is the hrst Lyapunov exponent. 



(H3). A(i? e (A(/.)))| M=WJ =d^0 



In this case, checking condition (HI) is not as easy as condition (SN1) 



Different methods are available in order to compute Hopf bifurcation points 
(see HI])- Our algorithm is based on computing the bialternate product 
of the Jacobian matrix of the system 2 J(j, P) Id where Id is the identity 
matrix. It is known that the determinant of this matrix vanishes if and only 
if the Jacobian matrix has two opposed eigenvalues (its eigenvalues are A, + Xj 
where A; and Xj are eigenvalues of the initial matrix J(j,P)). Therefore at 
these points where the bialternate product vanishes we have to check that there 
exists a purely imaginary eigenvalue to avoid cases where the system has two 
real opposed eigenvalues. Even with this step, the bialternate product method is 
much more efficient than other methods based on the characteristic polynomial 
such as Kubicck's method l32l. 



A. 5 Bautin bifurcation 

If along the line of the Hopf bifurcations the first Lyapunov coefficient van- 
ishes and a subcritical Hopf bifurcation becomes supercritical when varying the 
parameters, under the differential conditions of [331 theorem 8.2], the system un- 
dergoes a Bautin bifurcation. At the points where the first Lyapunov coefficient 
vanishes, the algorithm numerically checks these conditions. 
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